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INTRODUCTION 

This  project  focuses  on  a  general  problem  encountered  in  the  application  of  low  coherence 
interferometry  (LCI,  also  referred  to  as  Optical  Coherence  Tomography,  OCT]  to  measure  the 
fundamental  optical  properties  of  scattering  materials,  in  the  present  case  of  porous  ceramic  TBCs 
and  other  materials  used  in  extreme  environments:  namely,  that  currently  available  theoretical 
models  describe  experimental  data  inadequately  and  that  better  models  are  needed.  This  problem  is 
in  part  application-driven,  namely  based  on  the  need  to  be  able  to  extract  the  radiative  properties 
from  the  shape  the  LCI  signal.  On  the  other  hand,  the  solution  to  this  problem  also  has  a  broader 
value,  since  it  enables  more  accurate  computation  of  the  radiative  heat  transfer  through  TBCs  and 
other  materials  used  in  extreme  environments. 

Objectives 

1.  The  main  objective  is  to  develop  an  analytical  model  of  radiative  transfer  through  porous 
ceramic  thermal  barrier  coatings.  The  model  to  be  developed  should  compute  the  integrated 
optical  path-length  as  a  function  of  several  basic  parameters  including  bulk  scatter  and 
absorption,  and  takes  into  account  a  multilayer  structure  and  interface  reflections. 

2.  The  second  objective  is  to  validate  the  analytical  model  with  random  number  simulations  under 
a  variety  of  conditions,  including  coating  thickness,  scatter  and  absorption  parameter  values  and 
interface  effects. 

3.  The  third  objective  is  to  demonstrate  that  relevant  TBC  parameters  can  be  measured  with  LCI, 
by  applying  the  model  to  experimental  LCI  data  obtained  on  a  range  of  test  materials. 

Change  in  Objectives:  none 

Status  of  Effort 

As  of  the  date  of  this  report, 

1.  A  refined  random  walk  simulation  code  (in  two  languages]  has  been  finished  and  tested 

2.  A  new  model  has  been  derived  with  which  the  integrated  pathlength  can  be  computed  of 
photons  traveling  through  an  isotropically  scattering  medium. 


3.  This  new  model  was  tested  with  the  random  walk  model. 

As  of  to  date  the  project  completion  is  as  shown  in  the  following  Table. 

Table  1:  project  status  summary 


Task 

Description 

Completion 

(%) 

Anticipated 

completion 

1 

Review  the  currently  available  theoretical  models 

100 

N/A 

2 

Develop  a  new  radiative  transport  model 

75 

June  2017 

3 

Apply  random  walk  model  to  test  model  approaches 

75 

June  2017 

4 

Apply  the  theory  to  experimental  data  on  TBCs 

20 

June  2017 

5 

Report  on  results  and  future  recommendations 

50 

July  2017 

Project  completion  is  considered  to  be  overall  well  over  50%  from  a  technical  point  of  view.  At  least 
400  hours  have  been  dedicated  to  the  project,  to  which  a  total  of  400  hours  were  assigned  in  the 
contract.  For  this  period,  280  hours  will  be  billed  to  the  contract. 
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DETAILED  DESCRIPTION  OF  WORK  PERFORMED  DURING  THIS  REPORT  PERIOD 
Task  1.  Review  the  currently  available  theoretical  models 


The  OCT  signal  typically  decays  as  a  function  of  L,  the  coherence  path-matching  separation  between 
signal  and  reference  beam  photons.  This  difference  in  path  length  (traveled  by  signal  and  reference 
beam  photons]  is  the  result  of  signal  photons  travelling  a  distance  L  through  the  sample.  Since  this 
distance  L  is  related  to  scatter  and  absorption  properties  of  the  material,  the  measured  decay 
function  contains  precise  information  on  the  optical  properties  of  the  material  under  inspection, 
which  is  the  purpose  of  OCT  measurement.  The  available  models  to  date  do  not  with  sufficient 
accuracy  describe  measured  OCT  signals  of  even  relatively  simple  systems.  For  more  complex 
systems,  such  as  TBCs  that  are  partially  infiltrated  with  CMAS,  current  models  fall  short  of  being 
able  to  extract  useful  data.  To  account  for  this  deficiency,  the  purpose  of  this  project  is  to  improve 
models  for  use  in  OCT  measurements  on  TBCs. 

At  the  start  of  the  project,  the  main  model  approaches  identified  were 

1.  Diffusion  theory  for  highly  scattering  media,  see  ref.  [1] 

2.  Makino  4-flux  type  models,  see  ref.  [2] 

3.  Kubelka-Munk  flux  model,  see  ref.  [3] 

4.  Scallan  flux  model,  see  ref.  [4] 

5.  Fresnel-Huygens  derived  theory,  ref.  [5] 

Two  of  the  models,  #1  and  #5,  specifically  contain  a  parameter  that  is  needed  in  the  LCI 
measurements,  namely  the  integrated  pathlength,  while  the  other  approaches  do  not.  In  this  Task,  it 
was  investigated  how  well  these  two  models  describe  simulation  and  real  data,  and  which 
modifications  might  be  needed  to  further  improve  them.  This  comparison  between  analytical 
models  was  made  for  a  significantly  simplified  system,  being  that  of  a  semi-infinite,  isotropically  and 
purely  scattering  medium;  i.e.  effects  due  to  absorption,  interface  reflections,  anisotropy,  and 
multilayer  were  not  included.  The  underlying  thought  is  that  if  the  model  cannot  describe  this 
simplified  case  well,  it  may  not  be  adequate  for  more  complex  structures  either. 


First,  the  diffusion  theory  by  Popescu  and  Dogariu1  has  as  the  central  equation  describing  the 
energy  flux  of  the  diffusively  scattered  photons  with  integrated  pathlength  L  (in  their  publication 
called  s]  as 


J (L)  =  (4nlt/3)  3/2z0vL  5/2exp 


(1) 


where  lt  is  the  mean  free  path,  z0  is  a  so-called  extrapolation  length  that  carries  information  on  the 
effective  reflectivity  at  the  boundaries,  and  v  is  the  average  transport  velocity.  It  is  noted  that /(L) is  not 
the  same  as  the  probability  density  P(L)  of  optical  path  lengths  through  the  medium.  Unfortunately, 
upon  consultation  of  various  publications  by  the  authors,  it  is  not  entirely  clear  how  the  above 
equation  was  derived,  nor  how /(L)  relates  to  P(L).  It  appears  that  the  above  theory  is  obtained  by 
work  from  Patterson  et  al,6  and  from  various  textbooks.  Furthermore,  it  is  claimed  that  the  above 
equation  describes  measured  data  well,  however  we  found  that  the  equation  is  in  fact  not  entirely 
adequate  for  both  experimental  and  simulation  data.  For  brevity,  the  above  model  is  not  discussed 
in  detail. 
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The  second  model  that  was  investigated  is  the  model  derived  by  Thrane  et  al  from  Fresnel-Huygens 
diffraction  theory.5  The  Thrane  model  defines  the  normalized  signal  current  as  a  function  of 
integrated  pathlength  L  (in  their  publication  reffered  to  as  z)  as 


(i2(Q) 
( i2)o 


=  exp(— 2  /rsL)  + 


4exp(-^sL)[l-exp(-//sL)] 

l+wl/wli 


+  [1  -  exp (— /isL)] 


2  ?ji 

2  f 


(2) 


where  /rs  is  the  scatter  coefficient  (in  units  inverse  length]  and  ws  and  wH  are  optical  parameters. 
Analysis  of  the  Trane  model  indicates  that  it  does  not  account  sufficiently  for  the  strongly  non¬ 
exponential  dependence  of  signal  with  integrated  pathlength,  found  both  in  simulation  and 
experimental  data.  It  appears  that  while  the  model  was  derived  to  account  for  both  single  and 
multiple  scattering,  it  is  more  strongly  favouring  single  scattering  than  is  witnessed  from  data. 


Originally,  and  as  outlined  in  the  proposal,  the  analysis  of  these  currently  available  models  was 
considered  to  be  sufficient  to  lead  to  an  improved  model.  The  original  plan  was  to  select  a  currently 
available  model  and  find  the  modifications  needed  to  better  describe  both  simulated  and  measured 
data. 


Since  thorough  analysis  of  the  above  models  did  not  lead  to  either  model  describing  simulated  and 
experimental  data  adequately,  and  it  was  not  clear  which  improvements  needed  to  be  made  to  make 
theory  and  data  agree  better,  it  was  decided  to  derive  a  new  model  from  first  principles.  This  effort 
was  in  part  motivated  by  the  fact  that  multiple  scattering  should  be  well  described  by  Poisson 
statistics,  and  that  such  statistical  processes  can  be  well  tested  with  random  number  simulations. 

It  is  noted  that  therefore  the  original  Task  2  description,  Modify  currently  available  models,  has 
been  replaced  by  Develop  a  new  optical  pathlength  model. 

Task  2.  Develop  a  new  optical  pathlength  model  and  Task  3.  Apply  random  walk  model  to  test 
model  approaches 

In  this  section,  we  will  describe  work  performed  in  both  tasks  together  for  convenience.  First,  we 
briefly  discuss  the  random  walk  model,  since  it  is  relevant  in  the  discussion  on  both  the 
development  of  the  new  model  (Task  2]  and  in  testing  it  (Task  3],  Some  model  parameter 
definitions  are  shown  in  Figure  1.  The  scattering  medium  is  an  infinite  slab  with  thickness  D. 
Photons  are  scattered  randomly  multiple  times  and  travel  a  total  (integrated]  pathlength  within  the 
medium  of  L  =2  Incoming  photons  are  considered  to  be  traveling  as  a  pencil  ray,  as  shown, 
whereas  outgoing  photons  in  all  directions  are  collected. 


Figure  1.  Schematic  of  multiple  scatter  events  inside  a  TBC 


4 


Two  random  walk  simulation  codes  were  adapted  and  improved  upon,  one  written  in  Labview  and 
one  in  C++.  In  particular,  the  effects  of  integrated  path-length  needed  to  be  incorporated  in  the 
simulation  code,  and  tested.  The  validity  of  the  two  codes  was  tested  in  several  ways,  including  the 
predictions  of  total  transmission  and  reflection.  The  two  codes  were  also  compared  with  each  other, 
in  order  to  check  consistency.  For  instance,  Figure  2  shows  two  integrated  pathlength  distributions 
computed  with  both  models,  using  ^ s  =  200  mm1. 


Integrated  pathlength,  L(mm] 

Figure  2.  Comparison  of  pathlength  distributions,  as  computed  by  two 
random  number  simulation  codes  developed  for  this  project. 

Ideal  case:  infinitely  thick  slab  with  scatter  only 

First,  we  look  only  at  scattering  process  (i.e.,  absorption  is  set  to  zero],  and  furthermore  neglect 
reflections  at  the  interface.  All  diffusely  reflected  photons  are  considered,  i.e.  the  collection  optics 
are  assumed  to  be  irrelevant. 

The  first  questions  that  arises,  is  what  defines  an  infinite  slab  thickness,  and  what  does  the 
simulated  signal  look  like  in  such  a  case.  For  that  purpose,  scatter  coefficient  (/rs]  and  layer 
thickness  (Z)]  was  varied  in  a  set  of  simulations,  of  which  Figure  3  shows  an  example.  It  can  be  seen 
that  the  difference  in  D  manifests  itself  only  in  the  tail  pathlength  distribution.  For  instance,  the 
simulated  curves  for  (D=0.3  mm,  /rs= 20  mm1]  and  (D=0.6  mm,  /j.s= 20  mm1]  are  overlapping  up  till 
approximately  L  =  0.1  mm,  then  diverge  at  higher  values.  The  same  is  the  case  for  the  lines  for 
(D=0.3,  jirs=40]  and  (D=0.6,  jis= 40],  where  the  divergence  occurs  at  slightly  higher  values  of 
approximately  L  =  0.2  mm.  This  can  be  understood  to  be  the  result  of  the  photons  on  average 
traveling  a  longer  pathlength  in  a  thicker  medium  ,  whereas  for  short  integrated  pathlengths  the 
layer  thickness  is  irrelevant. 
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Integrated  pathlength,  L  (mm) 

Figure  3.  Effect  of  scatter  coefficient  and  layer  thickness  on  the  heterodyne  efficiency 

Another  set  of  integrated  pathlength  distribution  simulations  is  shown  in  Figure  4,  where  the  slab 
thickness  D  was  varied  in  steps  of  factor  2  while  the  scatter  coefficient  /j.s=  200  mm1.  Photons  that 
reach  the  boundary  were  allowed  to  escape  at  the  other  end.  Again,  the  effect  of  D  on  the  pathlength 
distributions  shows  up  in  the  tail  of  L,  while  for  small  values  the  curves  overlap. 


L  (mm) 


Figure  4.  Integrated  pathlength  distributions  for  different  values  of  slab  thickness  D 

On  the  other  hand,  the  difference  in  /rs  values  manifests  itself  predominantly  in  the  lower  values  for 
for  L  in  the  pathlength  distribution,  as  is  shown  in  Figure  5. 
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L[mm] 

Figure  5.  Integrated  pathlength  distributions  as  a  function  of  scatter  coefficient  /rs. 

Thus,  an  infinitely  thick  coating  exhibits  a  near  straight  line  in  the  tail  on  a  log-log  plot,  which 
indicates  a  power-law  functionality.  From  these  simulations,  a  condition  for  what  defines  an 
infinitely  thick  coating  is  therefore  approximately  D  x  >  100.  The  slope  of  the  tail  in  the  log-log 
plot  is  not  affected  by  /rs,  only  the  low  L  value  end  of  the  pathlength  distribution. 

New  optical  pathlength  distribution  model 

In  this  section,  the  basics  of  a  new  analytical  model  derived  as  part  of  the  AFOSR  project  will  be 
discussed.  In  parallel,  OCT  model  work  has  also  been  undertaken  at  MetroLaser,  Inc.  (Irvine, 
California],  under  a  USAF  SBIR  contract.  In  brief,  the  approach  taken  by  MetroLaser  is  as  follows.  It 
was  found  that  the  OCT  signal  may  be  approximated  by  (and  fit  with]  a  power  law  function  with 
form 

/l(L)oc(l+i)"”,  (3) 

where  Ls  is  a  characteristic  length  scaling  factor  and  n  expresses  the  degree  of  non-exponentiality 
and  is  considered  to  relate  to  the  number  of  scatter  events  contributing  to  the  signal.  This  decaying 
signal,  in  fact  any  decaying  signal,  may  be  defined  as  a  sum  of  exponential  components.  In  the 
continuous  case,  this  is  defined  by  the  Laplace  transform 

A  (L)  =  J“  D  (x)  exp  (-^)  dx,  (4] 

Where  D(x)  is  the  amplitude  distribution  function  (or  probability  density  function]  for  the  various 
values  of  some  characteristic  pathlength  x.  When  performing  an  inverse  Laplace  transform  on  this 
function,  it  was  demonstrated  by  MetroLaser  scientists  that  this  results  in  a  Poisson-like  function 

D(x)0C&)n  leXP(“x)’  W 

which  may  be  used  to  correlate  the  signal  with  scatter  coefficients  of  the  medium. 

From  the  point  of  view  of  signal  processing,  this  approach  is  robust  since  the  inverse  Laplace 
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Transform  is  performed  on  a  noiseless  dataset,  i.e.,  the  fit  with  Eq.  (3],  That  the  kernel  is  a  Poisson 
distribution  appears  in  accordance  with  diffusion  of  light  through  a  scattering  medium,  although  the 
physical  meaning  of  the  exponent  n  is  not  entirely  clear  yet. 

As  mentioned  above,  the  purpose  of  the  work  discussed  in  the  present  report  is  to  likewise  provide 
further  theoretical  underpinning  of  OCT  measurements.  The  approach  is  to  derive  equations  from 
first  principles  and  compare  the  theoretical  results  with  OCT  measurements  and  numerical 
simulations.  It  is  therefore  of  interest  as  well  how  the  resulting  theory  compares  with  the  above 
approach.  In  addition,  OCT  measurements  at  MetroLaser  as  a  function  of  concentration  [C]  of 
scatterers  indicate  that  the  peaks  of  the  distribution  functions  scale  by  [C]-1/3,  apparently  violating 
Beer's  law  for  optical  transmission,  while  the  random  walk  simulations  predict  a  scaling  with  [C]1, 
in  accordance  with  Beer's  law.  One  way  to  find  out  how  this  difference  may  arise  is  to  derive  a 
theoretical  framework  and  test  various  aspects  of  it  with  random  walk  models.  Finally,  the 
theoretical  framework  should  allow  us  to  evaluate  under  what  conditions  certain  approximations 
are  valid. 

In  deriving  an  analytical  set  of  equations  for  the  integrated  pathlength,  the  approach  taken  in  the 
present  work  was  to  define  the  problem  in  terms  of  probability  distribution  functions  for  the 
number  of  photons  that  travel  a  total  pathlength  L  through  a  homogenously  scattering  medium  with 
scatter  coefficient  /rs  following  n  scatter  events.  For  instance,  one  can  evaluate  the  first  passage 
problem,  i.e.,  evaluate  the  number  of  steps  required  to  escape  once  the  particle  has  entered  the 
medium.  The  problem  can  be  sub-divided  into  deriving  probability  functions  for  the  number  of 
photons  that  travel  a  distance  L  while  undergoing  n  =  1,  2,  3,  ...  scatter  events.  The  summation  of 
these  probability  functions  should  correspond  to  the  OCT  signal  beam;  i.e. 

PtotaliQ  =  -Pi  00  +  P2(L)  +  -  =  EnPn(P)-  (6) 

It  was  found,  by  simulation,  that  the  integrated  individual  terms  over  the  path-length,  i.e., 

Pn  =  fPn(L)dL.  (7) 

were  independent  of  scatter  cross-section  jxs,  as  shown  in  Figure  6.  If  true,  this  somewhat 
surprising  result  should  also  follow  from  deriving  the  explicit  equations  for  Pn(L).  The  values  for  n 
up  to  10  are  shown  in  the  table  to  the  right. 


n 

Pn 

1 

0.182 

2 

0.109 

3 

0.074 

4 

0.055 

5 

0.043 

6 

0.035 

7 

0.029 

8 

0.024 

9 

0.021 

10 

0.018 

Figure  6.  Probabilities  P„,  from  random  number  simulations. 

In  the  following  paragraphs,  the  derivation  of  these  terms  Pn(L)  will  be  discussed,  starting  with  a 
derivation  of  the  probability  distribution  function  Pi(L)  and  hence  Pt. 
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PiXL)  derivation 

Consider  the  following  geometrical  arrangement,  see  Figure  7,  whereby  a  collinear  low  coherence 
beam  (shown  as  a  red  arrow]  is  considered  to  be  incident  onto  the  sample  at  normal  incidence  to 
the  surface  (shown  as  a  circle]. 


Figure  7.  Optical  geometry  of  one  particular  n  =  1  scatter  event. 

The  photons  travel  into  the  sample  along  the  same  direction,  and  scatter  at  some  position  p1  along 
this  line,  at  a  distance  lt  from  the  boundary.  The  location  is  not  a  fixed  position:  the  probability 
density  function  of  the  number  of  photons  scattered  at  a  position  is  defined  by  the  exponential 
function 

Ppi(li)  =  PsexpC-pA)  .  (8] 

Upon  random  scattering  at  location  p1;  the  photon  travels  in  some  arbitrary  direction.  Suppose  the 
photon  travels  in  the  direction  indicated  by  the  arrow  in  Figure  7,  is  reach  the  boundary  and 
escapes,  and  hence  is  scattered  only  once.  The  distance  to  the  boundary  is  defined  as  lb,  such  that 
the  total  distance  that  the  photon  has  traveled  is  L  =  +  lb.  Then,  the  probability  distribution  of 

total  distance  travelled  (L]  for  the  case  of  singly  scattered  photons,  Pi(L),  may  be  defined  as  the 
integral  product 

PliD  =  C  Pescdb  =  L  -  IDPpiOOdh  ,  (9] 

where  Pesc(Zb  =  L  —  h)  is  the  escape  probability  for  a  photon  traveling  a  distance  lb  from  the  scatter 
location  p1  to  the  boundary,  i.e.,  under  the  condition  that  it  has  traveled  a  path  with  total  length  L. 

The  following  equation  was  derived  for  the  escape  probability  as  a  function  of  distance  db  to  the 
boundary: 

P»c(4)  =  (/lexp(-^)d8.  (10) 

In  the  above  example  if  single  scattering  photons,  db  is  simply  equal  to  A  closed  form  for  the 
above  integral  may  be  obtained,  but  is  not  very  compact  and  hence  not  shown  here  for  brevity.  Also, 
note  that  the  integral  has  a  singularity  for  8  =  -,  nevertheless  Pesc(db )  goes  to  zero  as  8  goes  to 
The  analytical  solution  of  Eq.  (10]  and  random  number  simulation  results  are  compared  in  Figure  8. 
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djj  (mm) 


Figure  8.  Escape  probability  Pesc(db ),  computed  with  random  number  simulations 

and  with  Eq.  (10). 

When  combining  Eqs.  (8)-(10),  one  can  derive  the  following  equation  for  P-^L): 

pi (L)  =  (i  _  MsexP (-/isL)  ■  (11) 


The  total  number  of  photons  that  are  scattered  only  once  can  thus  be  found  to  be 

P1  =  SP1(.L)dL  =  (j-±)~  0.182. 


(12) 


It  is  thus  found  that  P1  is  a  constant  and  does  not  depend  on  the  scatter  coefficient.  The  value  of  0.18 
was  verified  with  simulations,  shown  in  the  table  in  Figure  6. 

higher  order  Pn(L)  terms 

The  first  step  in  the  derivation  of  the  higher  order  terms  is  to  define  the  probabilities  of  traveling  a 
total  distance  L  while  undergoing  2,  3,  etc.  scatter  events.  This  probability  can  be  defined  as  the 
integral  over  three  probabilities:  the  probability  of  scatter  at  location  p1;  the  probability  of  scatter  at 
p2  and  the  escape  probability  from  p2 under  the  condition  that  L  =  lt  +  l2  +  lb: 

p2 a)  =  f  f  PplGl)Pp2(h)Pesc(db)dPldp2  .  (13) 


The  geometrical  arrangement  of  one  such  n  =  2  scatter  event  is  shown  in  Figure  9. 


Figure  9.  Optical  geometry  of  one  particular  n  =  2  scatter  event. 
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The  solution  to  Eq.  (13]  is  non-trivial,  and  full  derivation  will  be  presented  at  a  later  stage.  Part  of 
the  solution  is  the  definition  of  the  probability  of  traveling  a  distance  dz  along  the  z-axis  following 
an  isotropic  scatter  event,  which  was  derived  to  be  described  by: 


p  (d  }  rvexP(-^sdzsin  10) 

dz(~  z '  2n'o  sin0 


d6 


(14] 


The  following  graph  shows  the  probability  density  function  of  PdZ(dz)  together  with  a  random 
number  simulation  of  the  same  process: 


Figure  10.  probability  of  traveling  a  distance  dz  following  a  scatter 
event:  random  number  simulation  versus  Eq.  (14]. 


Working  model  of  higher  order  P„(L)  terms 

While  explicit  forms  of  Pn(L)  are  preferred,  the  following  semi-theoretical  approach  was  used  to 
derive  a  working  model. 

In  general,  the  probability  of  undergoing  n  scatter  events  over  a  distance  L,  with  an  average  number 
of  events  A  =  p.sL,  is  given  by  the  Poisson  distribution 

Ps(n,  L)  =  C^expW)  (15] 

The  random  walk  simulation  was  used  to  test  how  this  distribution  comes  into  play  in  the  overall 
pathlength  distribution  signal  Ptota;(L ).  It  became  evident  that  instead  of  the  above  equation,  a 
slightly  different  version  was  fully  describing  the  number  of  photons  that  travel  a  distance  z,  let's 
call  it  F  : 

Fn(L)  =  constant  x  ,  (16] 


where  the  constant  depends  on  n.  Specifically,  the  function  that  was  used  to  fit  the  histogram  curves 
from  random  number  simulations  was 


Pri(P)  Q-nF-s 


insL)n  1exp(-^sL) 
(n-1)! 


(17] 
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where  an  are  scaling  factors.  Figure  11  shows  four  histograms  for  different  values  of  n,  using 
/4=100  mm1  (the  actual  value  of  fis is  not  relevant],  together  with  the  fitting  function  of  Eq.  [17], 
Note  that  n.  L.  and  /rs  in  fitting  function  Fn (D  are  constants:  only  the  scaling  factor  an  is  used  as  a 
fitting  variable. 

Two  observations  are  made;  firstly,  for  n  =  1  Eq.  [17]  reduces  to  Fi(L)  =  anjusexp(— nsL)  and  the 
fitting  value  %  is  found  to  be  0.18,  i.e.  close  to  the  value  of  Q  —  of  Eq.  (12]. 

Secondly,  the  other  fitting  parameters  a3,  a5  and  a10  as  shown  in  Figure  11  are  very  close  to  the 
values  P3,  P5  and  Pio  as  shown  in  the  table  of  Figure  6.  Observation  of  other  results  revealed  that  in 
general  an  has  approximately  the  same  value  as  Pn.  One  can  see  that  Eq.  (17]  divided  by  an  are 
probability  functions  with  unity  total  probability,  hence  an  are  simply  the  normalization  constants 
that  determine  the  integrated  values  of  the  probability  density  functions. 


Figure  11.  Probability  distribution  functions  P„(L]  from  random  number  simulation 
and  using  fitting  function  Eq.  (17],  for  n  =  1,  3,  5  and  10.  In  all  cases,  ^s=100  mm1 

Thus  one  can  assume  that  the  fitting  function  Eq.  (17]  serves  as  a  good  working  model  for  the  OCT 
signal  modeled  by  the  random  number  simulation  code.  The  extent  to  which  this  model  is  suitable 
for  real  OCT  signals  will  be  evaluated  subsequently.  So  the  signal  beam  photon  count  in  this 
idealized  case  (since  the  optical  geometry  is  not  taken  into  account,  and  several  other  idealizations 
are  in  effect]  can  be  approximated  by  rewriting  the  terms  of  Eq.  (17]  and  summing  over  n : 

Signal  oc  S(L)  =  Z“=1Fn(L)  =  Zn=i^^n_1exp(-/isL)  (18] 
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The  upper  limit  in  the  summation  is  taken  to  be  infinite,  however  in  reality  there  is  likely  an  upper 
limit  in  the  number  of  scatter  events  that  can  be  measured.  However,  the  decrease  in  coefficients  an 
already  take  care  of  this  reduced  sensitivity  of  the  measurement  to  higher  numbers  of  scatter 
events.  Since  the  exponential  term  does  not  contain  n,  and  to  start  the  summation  with  n  =  0,  we  can 
rewrite  Eq.  (18]  as 

S(L)  =  Msexp(-/rsL)£“=0^(/rsL)n  (19] 

An  example  of  a  comparison  between  Eq.  [19]  and  the  full  histogram  as  computed  with  RN 
simulation  is  shown  in  Figure  12,  which  was  obtained  for  /rs=100  mm1.  Instead  of  an  infinite 
number  of  terms,  the  graph  shows  the  summation  of  Eq.  [19]  over  terms  with  n  =  0  -  130.  As  might 
be  expected,  the  power  series  of  Eq.  [19]  was  found  to  describe  the  simulated  signal  histogram  well, 
for  a  given  range  of  L:  the  smaller  the  number  of  terms  (i.e.,  for  smaller  upper  n),  the  shorter  the 
range  of  L  over  which  Eq.  [19]  compares  well  with  the  RN  simulation.  The  range  of  L  which  is  well 
described  by  Eq.  (19]  is  roughly  nu Pper/ps- 


Figure  12.  Signal  function  as  obtained  from  RN  calculations  and  computed  with  Eq. 

(19]  for  a  finite  range  of  terms,  for  n  =  0  to  130,  for  /is=100  mm1. 

It  is  noted  that  the  signal  function  has  the  form  of  a  Poisson  Generating  Function,  i.e. 

S(L)  =  PG(an;x)  =  (20] 

The  use  of  an+1  instead  of  an  in  Eq.  (19]  is  arbitrary  and  hence  insignificant. 


Relation  with  MetroLaser  model 

In  this  section,  the  goal  is  to  see  how  the  above  derivation  relates  to  the  MetroLaser  (ML]  model. 

The  basic  principles  of  the  ML  model  is  that  the  OCT  signal  of  a  randomly  (diffuse]  scattering 
medium  appears  to  fit  well  to  a  power  law  function  of  the  form 


Ppl (£) 


a 

(1 +bL)c 


(21] 


In  turn,  by  means  of  the  inverse  Laplace  Transform  this  power  law  function  can  be  converted  to  a 
Poisson-like  kernel  function,  from  which  a  correlation  with  scatter  coefficient  may  be  obtained. 
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First,  it  is  worth  checking  how  well  the  simulated  signal  is  fit  with  Eq.  (21].  This  can  be  seen  in 
Figure  13.  As  can  be  seen,  Eq.  [21]  doesn't  fit  the  simulated  data  entirely  well.  The  fit  values 
obtained  were  a  =  1.01;  b  =  24.0  and  c  =  1.82.  In  order  to  make  a  better  fit,  a  baseline  value  can  be 
added,  as  is  done  in  the  ML  approach.  However,  the  physical  interpretation  of  this  baseline  value  is 
not  clear  in  the  present  model,  since  it  predicts  that  the  signal  decreases  to  zero  for  large  enough  L. 


Figure  13.  Signal  function  S(L]  and  a  fit  with  power  law  function  from  Eq.  (21],  for 
/rs=100  mm1. 

Next,  instead  of  starting  with  a  signal  (either  simulated  or  real]  and  trying  to  extract  a  kernel 
function  via  the  inverse  Laplace  Transform,  it  is  of  interest  to  see  how  the  signal  function  Eq.  (19] 
itself  can  be  expressed  as  a  (direct]  Laplace  Transform,  and  if  so,  to  what  function  this 
transformation  leads. 

In  general,  any  power  series 

2n=o  bnxn  (22] 


can  be  converted  to  an  integral  (i.e.,  continuous]  form  of  the  Laplace  Transform  by  making  three 
substitutions:  (1]  exchange  integer  n  for  a  continuous  parameter  t;  (2]  defining  — s  =  lnx;  and  (3] 
multiplying  bn  with  the  Dirac  delta  function  S(t  —  ri).  In  that  case,  one  gets 

£“=o  bnxn  =  J”  h(t)exp(— st)dt  =  £{6(t)}.  (23] 


Thus,  we  see  that  the  signal  function  S(L]  can  be  written  as  a  Laplace  Transform  of  a  function  that 
contains  the  constants  an+1: 

S(L )  =  £{a(t)}  (24] 

where 

a(t)  =  {ts^rexp(-/isL)  S(t-n)  (25] 


Thus,  we  see  that  in  the  current  approach,  the  Poisson  generating  signal  function  Eq.  (19]  can  in 
principle  be  converted  into  a  Laplace  transform.  Although  numerical  values  of  the  coefficients  an  are 
known,  it  will  be  of  value  to  define  them  analytically,  which  is  being  investigated.  Whether  this  leads 
exactly  to  a  power  law  function  of  the  form  of  Eq.  (21]  remains  to  be  seen.  Considering  that  Eq.  (21] 
does  not  fit  simulated  data  perfectly  well  unless  a  baseline  is  added,  the  current  approach  may  likely 
lead  to  a  different  signal  function. 
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CHANGES  IN  SCHEDULE 

A  no-cost  extension  of  will  be  requested  for  this  project,  with  a  new  project  end  date  of  30  June 
2017. 

WORK  PLAN  FOR  NEXT  REPORTING  PERIOD 

In  the  next  reporting  period,  the  remaining  tasks  of  this  project  will  be  conducted,  including 

1.  Finalize  the  new  analytical  model  by  incorporating  the  effect  of  absorption  and  interface 
reflections  on  the  integrated  pathlength  distribution 

2.  Evaluate  the  new  model  with  random  walk  simulations 

3.  Analyse  experimental  data  with  the  new  model 

4.  Submit  the  final  report. 

PERSONNEL  SUPPORTED 

Personnel  working  on  this  project  was  B.  Heeg,  Ph.D. 

COST  SUMMARY 

As  of  this  date,  280  hours  have  been  billed  to  the  project,  leaving  a  total  of  120  hours  to  bill  to  the 
project. 
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